/*
Author: Kuljit Virk
Date  : 22 Jan 2012
Purpose: Input output for the superlattice code

Read in: Parameter file for tight-binding Hamiltonians
*/

#include <iostream>
#include <string.h>
#include <complex>
#include <math.h>
#include "mkl_types.h"
#define MKL_Complex16 std::complex<double>
#include "mkl.h"
#include "mkl_lapacke.h"

typedef MKL_Complex16 complex16;
typedef MKL_INT integer;

using namespace std;


#define INTERFACE 1
#define BULK 2

struct tb{
  float Esa, Epa, Esc, Epc, Essa, Essc;
  float Vss, Vxx, Vyy,Vxy,Vscpa, Vsapc, Vssapc, Vpassc;
  float Usoa,Usoc;
};

enum material {GaAs, AlAs};



void get_tb_parameters_vogl(struct tb params[],const int so){
/* For the conventional tight binding model, we use
   10-orbital model for each type of atom. We only
   address zinc-blende crystals, and the form of the
   Hamiltonian is taken from Vogl, and we use
   Chadi, and Chadi and Cohen papers to specify the
   spin-orbital part of the Hamiltonian.

   Parameters are:
   E(s ,a) E(p,a) E(s,c) E(P,c) E(s*,a) E(s*,c)

   V(s,s) V(X,X) V(X,Y) V(sa,Pc) V(sc,pa) V(s*a,Pc) V(Pa,s*c) 

   Here a,c label the anion and the cation respectively.
   In the case of GaAs, a = As, and c = Ga

   We arrange the orbitals differently from Vogl because
   we decompose our structure in layers. Each atomic plane
   is either 'a' or 'c' type. Within each the 10 orbitals are
   arranged as follows (numbers on the left enumerate them)

   For 'a':
   0:   |s  +>       
   1:   |px +>      
   2:   |py +>       
   3:   |pz +>      
   4:   |s* +> 
   5:   |s  ->       
   6:   |px ->      
   7:   |py ->       
   8:   |pz ->
   9:   |s* +> 

   For 'c':  numbers 10 .. 19

   Thus the matrix elements have the following enumeration

   E(s ,a) E(px..pz,a) E(4,a) 
   E(0 ,a) E(1....3,a) E(4,a)

   V(s,s) V(X,X) V(X,Y) V(sa,Pc)     V(sc,pa)   V(s*a,Pc) V(Pa,s*c) 

   Spin up
   V(0,10) V(1,11) V(1,2) V(0 ,11..13) V(10,1..3) V(4,11..13) V(1..3,14) 

   Spin down
   V(5,5) V(6,6) V(6,7) V(5 ,16..18) V(15,6..8) V(9,16..18) V(6..8,19) 

   SO coupling occurs within the cation or the anion
   V(px+,pz-)    V(py+,pz-)    V(px+,py+)
   Uso           iUso         -iUso 

   The coupling between atomic planes:  V_{As, Ga} = V(a,c). These
   components are:
   
   V(sa,Pc)   V(sc,pa)   V(s*a,Pc) V(Pa,s*c) 

   The file is specified in terms of columns, with one column per material
   and rows corresponding to the matrix element

*/
  struct tb p;
/* Vogl: GaAs
 */ 
  p.Esc = -8.3431;
  p.Epa =  1.0414; 
  p.Esa = -2.6569;
  p.Epc =  3.6686;
  p.Essa=  8.5914;
  p.Essc=  6.7386;
  
  p.Vss   = -6.4513;
  p.Vxx   =  1.9546;
  p.Vxy   =  5.0779;
  p.Vsapc =  4.3800;
  p.Vscpa =  5.7839;  // no
  p.Vssapc=  4.8422;  
  p.Vpassc=  4.8077;  // no

  // Taken from Vogl who quotes Chadi
  // The quoted values are splittings Da or Dc
  // The matrix element is splitting/3
  p.Usoa = 0.421/3.0 * so;
  p.Usoc = 0.174/3.0 * so;

  params[GaAs]=p;

/* Vogl: AlAs
 */ 
  p.Esa = -7.5273;
  p.Epa =  0.9833; 
  p.Esc = -1.1627;
  p.Epc =  3.5867;
  p.Essa=  7.4833;
  p.Essc=  6.7267;
  
  p.Vss   = -6.6642;
  p.Vxx   =  1.8780;
  p.Vxy   =  4.2919;

  p.Vsapc =  5.1106;

  p.Vscpa =  5.4965;

  p.Vssapc=  4.5216; 
  p.Vpassc=  4.9950; 

  // Taken from Vogl who quotes Chadi
  // The quoted values are splittings Da or Dc
  // The matrix element is splitting/3
  p.Usoa = 0.421/3.0 * so;
  p.Usoc = 0.024/3.0 * so;

  params[AlAs]=p;
  return;
};


void get_tb_parameters(struct tb params[],const int so){
  struct tb p;
/* Kimeck: GaAs: SO included
 */ 
  p.Esa = -3.53284;
  p.Epa =  0.2772 ;
  p.Esc = -8.11499;
  p.Epc =  4.57341;
  p.Essa=  12.33930;
  p.Essc=  4.31241;
  
  p.Vss   = -6.87653;
  p.Vxx   =  1.33572;
  p.Vxy   =  5.07596;
  p.Vsapc =  2.85929;
  p.Vscpa =  11.09774;  // no
  p.Vssapc=  6.31619  ;
  p.Vpassc=  5.02335;  // no
  p.Usoa  =  0.32703;
  p.Usoc  =  0.12000;

  params[GaAs]=p;

/* Klimeck: AlAs
 */ 
  p.Esa = -3.21537;
  p.Epa =  0.09711 ;
  p.Esc = -9.52462;
  p.Epc =  4.97139;
  p.Essa=  12.05550;
  p.Essc=  3.99445;
  
  p.Vss   = -8.84261;
  p.Vxx   = -0.01434;
  p.Vxy   =  4.25949;
  
  p.Vsapc =  2.42949;

  p.Vscpa =  13.20317;

  p.Vssapc=  5.83246;
  p.Vpassc=  4.60075;

  p.Usoa = 0.29145;
  p.Usoc = 0.03152;
  
  params[AlAs]=p;
  return;
};

inline complex16 operator *(double d, complex16 f){
  return complex16(real(f)*d,imag(f)*d);
};


void tb_Hamiltonian(complex16 *H, double k[], struct tb params){
  int i,j;
  int tbdim=20;
  int halfdim=10;

  for (int i=0;i<tbdim;i++)
    for (int j=0;j<tbdim;j++){
      H[i*tbdim + j]=0.0;
    }

  complex16 g[4];

  g[0] = cos(k[0]/4)*cos(k[1]/4);
  g[1] = complex16(0.0,sin(k[0]/4)*cos(k[1]/4));
  g[2] = complex16(0.0,cos(k[0]/4)*sin(k[1]/4));
  g[3] = -sin(k[0]/4)*sin(k[1]/4);

  cout << "g[0..3]"<<"\n";
  for (int i=0;i<4;i++)
    cout << g[i]<<"\n";

  // Diagonal terms
  H[0*tbdim+0] = params.Esa;
  H[1*tbdim+1] = params.Epa;
  H[2*tbdim+2] = params.Epa;
  H[3*tbdim+3] = params.Epa;
  H[4*tbdim+4] = params.Essa;

  // Diagonal terms
  H[10*tbdim+10] = params.Esc;
  H[11*tbdim+11] = params.Epc;
  H[12*tbdim+12] = params.Epc;
  H[13*tbdim+13] = params.Epc;
  H[14*tbdim+14] = params.Essc;

  // Off-diagonal terms
  H[0*tbdim+10] = complex16(params.Vss) * g[0];

  for (i=1;i<=3;i++){
    // V(sa, pc)
    H[0*tbdim+10+i] = complex16(params.Vsapc) * g[i];
    
    // V(pa, sc)
    H[i*tbdim+10]   = -complex16(params.Vscpa) * g[i];
  }

  // V(pa,pc)
  complex16 gg[3][3];
  gg[0][0] = g[0];
  gg[0][1] = g[3];
  gg[0][2] = g[2];

  gg[1][0] = g[3];
  gg[1][1] = g[0];
  gg[1][2] = g[1];

  gg[2][0] = g[2];
  gg[2][1] = g[1];
  gg[2][2] = g[0];

  for (i=1;i<=3;i++){
    for (j=1;j<=3;j++){
	// H[px_a,px_c] etc
	// H[px_a,py_c] etc
      if (i==j)
	H[i*tbdim+10+i] = params.Vxx * gg[i-1][j-1];
      else
	H[i*tbdim+10+j] = params.Vxy * gg[i-1][j-1];
    }
  }  
  // V(s* a,s* c) = 0 in Vogl paper
  H[4*tbdim+14] = 0.0;
  // V(s* a,p c) 
  H[4*tbdim+11] = params.Vssapc * g[1];
  H[4*tbdim+12] = params.Vssapc * g[2];
  H[4*tbdim+13] = params.Vssapc * g[3];

  // V(p a,s* c) 
  H[1*tbdim+14] = -params.Vpassc * g[1];
  H[2*tbdim+14] = -params.Vpassc * g[2];
  H[3*tbdim+14] = -params.Vpassc * g[3];

  // SO terms:
  H[1*tbdim + 3 + tbdim/4] =  params.Usoa * g[2];                 // < px+ | Hso | pz- >
  H[2*tbdim + 3 + tbdim/4] =  complex16(0.0,params.Usoa) * g[1];  // < py+ | Hso | pz- >
  H[1*tbdim + 2]           = -complex16(0.0,params.Usoa) * g[3];  // < px+ | Hso | py+ >

  H[11*tbdim + 13 + tbdim/4] =  params.Usoc * g[2];                // < px+ | Hso | pz- >
  H[12*tbdim + 13 + tbdim/4] =  complex16(0.0,params.Usoc) * g[1]; // < py+ | Hso | pz- >
  H[11*tbdim + 12]           = -complex16(0.0,params.Usoc) * g[3]; // < px+ | Hso | py+ >


  //----------- COPYING   
  // Spin down terms: anion
  for (i=0;i<tbdim/4;i++)
    for (j=i;j<tbdim*3/4;j++){
      H[(i+tbdim/4)*tbdim + j+tbdim/4] = conj(H[i*tbdim+j]);
    };

  // Spin down terms: cation
  for (i=tbdim/2;i<tbdim*3/4;i++)
    for (j=i;j<tbdim*3/4;j++){
      H[(i+tbdim/4)*tbdim + j+tbdim/4] = conj(H[i*tbdim+j]);
    };


  // Conjugate terms
  for (i=0;i<tbdim;i++)
    for (j=i;j<tbdim;j++){
      H[j*tbdim+i] = conj(H[i*tbdim+j]);
    }
};

void tm_template(complex16 *T, const int N){
  int N1 = N/2;
  int i,j;

  for (i=0;i<N;i++)
    for (j=0;j<N;j++){
      T[i*N+j] = 0.0;
    }
  
  for (i=0;i<N1;i++){
    T[i*N+i+N1]   = -1.0;
    T[(i+N1)*N+i] =  1.0;
  }
};

void layer_transfer_matrix(complex16 T[],const double E,complex16 H[], const int tbdim, const int layer_type){
  int dim1=tbdim/2;
  complex16 Tbb[tbdim*tbdim];
  complex16 Taa[tbdim*tbdim];
  
  complex16 Vab[dim1*dim1];
  complex16 Vba[dim1*dim1];
  complex16 Haa[dim1*dim1];
  complex16 Hbb[dim1*dim1];

  int i,j;
  for (i=0;i<dim1;i++){
    for (j=0;j<dim1;j++){

      Vab[i*dim1 + j]   = -H[i*tbdim + j + dim1];
      Vba[i*dim1 + j]   = -H[(i+dim1)*tbdim + j];

      Haa[i*dim1 + j]   = -H[i*tbdim + j];
      Hbb[i*dim1 + j]   = -H[(i+dim1)*tbdim + j + dim1];
    }
    Haa[i*dim1 + i ] += E;
    Hbb[i*dim1 + i ] += E;
  }

  // Set: Haa --> -inv(Vba)Haa using Lapack
  //      Solve:  Vba X = Haa  for X and set Haa ---> X  at the end
  MKL_INT info;
  MKL_INT ipiv[dim1];
  info = LAPACKE_zgesv(LAPACK_ROW_MAJOR, dim1, dim1, Vab, dim1, ipiv, Haa, dim1);
  if ( info != 0 ){
    cout << "zgesv for inv(Vab)*Haa failed, info = "<<info<<"\n";
    if (info > 0)
      cout << "U(i,i) for i = "<<info<<" is zero. NOT SOLUTION\n";
  }
  // Set: Hbb --> -inv(Vab)Hbb using Lapack
  //      Solve:  Vbb X = Hbb
  info = LAPACKE_zgesv(LAPACK_ROW_MAJOR, dim1, dim1, Vba, dim1, ipiv, Hbb, dim1);
  if ( info != 0 ){
    cout << "zgesv for inv(Vba)*Hbb failed, info = "<<info<<"\n";
    if (info > 0)
      cout << "U(i,i) for i = "<<info<<" is zero. NOT SOLUTION\n";
  }

  tm_template(Taa,tbdim);
  tm_template(Tbb,tbdim);

  for (i=0;i<dim1;i++){
    for (j=0;j<dim1;j++){
      Taa[i*tbdim+j] = -Haa[i*dim1+j];
      Tbb[i*tbdim+j] = -Hbb[i*dim1+j];
    }
  }
  
  // Replace Taa ---> Taa Tab
  complex16 alpha=1.0,beta=0.0;
  cblas_zgemm(CblasRowMajor,CblasNoTrans,CblasNoTrans,
  	      tbdim,tbdim,tbdim,&alpha,Taa,tbdim,Tbb,tbdim,&beta,T,tbdim);

  // Diagonalize
  /* complex16 eigval[tbdim], vl[tbdim*tbdim], vr[tbdim*tbdim];
  complex16 sum;
  if ( layer_type != INTERFACE ){
    info = LAPACKE_zgeev(LAPACK_ROW_MAJOR,'V','V',tbdim,T,tbdim,eigval,vl,tbdim,vr,tbdim);
    sum=0.0;
    for(int p=0;p<tbdim;p++){
      if (fabs(abs(eigval[p])-1.0) < 0.01)
	printf("%14.4e  %14.4e  %14.4e\n",-imag(eigval[p]),real(eigval[p]),abs(eigval[p]));
      sum += eigval[p];
    }
    cout << "sum = "<<sum<<"\n";    
  }
  */
  return;
};

void printeig(complex16 H[], const int& tbdim){
  int info;
  double energies[tbdim];
  info = LAPACKE_zheev(LAPACK_ROW_MAJOR,'V','U',tbdim,H,tbdim,energies);
  cout << "-----------------------------------------\n";
  for(int p=0;p<tbdim;p++)
    printf("%d  %14.4f\n",p,energies[p]);
  cout << "-----------------------------------------\n";
};
void printMatrix(complex16 H[], const int & tbdim);

int main(){
  cout << "Holy\n" ;
  struct tb params[2];
  complex16 *H;
  int tbdim = 20;
  double k[2];

  int info;
  complex16 T[tbdim*tbdim],vl[tbdim*tbdim],vr[tbdim*tbdim],eigval[tbdim];
  complex16 alpha=1.0,beta=0.0;
  double energies[tbdim];
  double pi = 4*atan(1.0);

  H = new complex16[tbdim*tbdim];

  k[0]=0.0;
  k[1] = 0.0;

  get_tb_parameters(params,0);
 
  tb_Hamiltonian(H,k,params[GaAs]);
 
  k[0]=0.0;
  cout <<"k = "<<k[0]<<endl;
  k[1]=0.0;

  get_tb_parameters(params,0);
  tb_Hamiltonian(H,k,params[GaAs]);
  printMatrix(H,tbdim);
  printeig(H,tbdim);
  cout <<"\n";

  tb_Hamiltonian(H,k,params[GaAs]);
  int n=tbdim/2;
  complex16 F[n*n];

  for (int i=0;i<tbdim/4;i++)
    for (int j=0;j<tbdim/4;j++){
      F[i*tbdim/2+j] = H[i*tbdim + j];
      F[(i+tbdim/4)*tbdim/2+j+tbdim/4] = H[(i+tbdim/2)*tbdim + j + +tbdim/2];
      F[(i)*tbdim/2+j+tbdim/4] = H[(i)*tbdim + j + +tbdim/2];
     }

  for (int i=0;i<tbdim/2;i++)
    for (int j=0;j<i;j++)
      F[i*tbdim/2+j]=conj(F[j*tbdim/2+i]);

  cout<<endl;
  printMatrix(F,tbdim/2);
  printeig(F,tbdim/2);
  printMatrix(F,tbdim/2);
  return 0;
  

  complex16 Tgaas[tbdim*tbdim];
  complex16 Talas[tbdim*tbdim];
  double E=1.0;

  for (int j=0;j<1;j++){
    k[0] = 0.50*j/4.0+0.1;
    k[1] = 1.0*j/4.0+0.1;
    tb_Hamiltonian(H,k,params[GaAs]);
    layer_transfer_matrix(Tgaas,E,H,tbdim,BULK);
    tb_Hamiltonian(H,k,params[AlAs]);
    layer_transfer_matrix(Talas,E,H,tbdim,BULK);
    // Replace Taa ---> Taa Tab
    cblas_zgemm(CblasRowMajor,CblasNoTrans,CblasNoTrans,
		tbdim,tbdim,tbdim,&alpha,Talas,tbdim,Tgaas,tbdim,&beta,T,tbdim);

    info = LAPACKE_zgeev(LAPACK_ROW_MAJOR,'V','V',tbdim,T,tbdim,eigval,vl,tbdim,vr,tbdim);
    for(int p=0;p<tbdim;p++)
      if (fabs(abs(eigval[p])-1.0) < 0.01)
	printf("%14.4e  %14.4e  %14.4e\n",-imag(eigval[p]),real(eigval[p]),abs(eigval[p]));
  }

  cout << "Done\n" ;
  
  return 0;
}


  /*
  cout << params[GaAs].Esa <<"\n";
  cout << params[GaAs].Epa <<"\n";
  cout << params[GaAs].Esa <<"\n";
  cout << params[GaAs].Epa <<"\n";
  cout << params[GaAs].Essa <<"\n";

  cout << params[AlAs].Esa <<"\n";
  cout << params[AlAs].Epa <<"\n";
  cout << params[AlAs].Esa <<"\n";
  cout << params[AlAs].Epa <<"\n";
  cout << params[AlAs].Essa <<"\n";
  */
void printMatrix(complex16 H[], const int & tbdim){
  printf("%1s"," ");
  for (int i=0;i<tbdim;i++){
    for (int j=0;j<tbdim;j++){

      if (abs(H[i*tbdim + j])>1.0e-4)
	if (real(H[i*tbdim + j]))
	  printf("%7.3f ",real(H[i*tbdim + j]));
	else
	  printf("%7.3fi",imag(H[i*tbdim + j]));
      else
	printf("%7s ","o");
      if (j==tbdim/4-1 || j==tbdim/2-1 || j==tbdim*3/4-1)
	printf("%1s","|");
    }
    cout << " \n ";
    if (i==tbdim/4-1 || i==tbdim/2-1 || i==tbdim*3/4-1){
      for (int ii=0;ii<tbdim;ii++)
	printf("%7s ","------");
      cout << " \n ";}      
  }
};
